P3704 [SDOI2017]数字表格
直接开始推式子,令 k=min(n,m)
ans=i=1∏nj=1∏mfgcd(i,j)=d=1∏kfdi=1∑nj=1∑m[gcd(i,j)=d]
上面那坨非常明显的莫反
=d=1∏kfdd∣T∑μ(dT)⌊Tn⌋⌊Tm⌋=T=1∏kd=1∏kfdμ(dT)⌊Tn⌋⌊Tm⌋=T=1∏k(d=1∏kfdμ(dT))⌊Tn⌋⌊Tm⌋
中间那坨可以直接 O(nlogn) 预处理,然后数论分块即可
代码
#include<iostream> #include<cstdio> #define MOD (1000000007) using namespace std; const int N = 1e6 + 5; typedef long long ll;
bool st[N]; int prime[N], mu[N], tot, n, m, g[N], f[N];
int qmi(int a, int b) { int res = 1; while (b) { if (b & 1) res = (ll)res * a % MOD; a = (ll)a * a % MOD; b >>= 1; } return res; }
void init() { g[0] = f[1] = mu[1] = g[1] = 1; for (int i = 2; i < N; i++) { g[i] = 1; f[i] = ((ll)f[i - 1] + f[i - 2]) % MOD; if (!st[i]) prime[++tot] = i, mu[i] = -1; for (int j = 1; i * prime[j] < N; j++) { st[i * prime[j]] = true; if (i % prime[j] == 0) break; mu[i * prime[j]] = -mu[i]; } } for (int i = 1; i < N; i++) { for (int j = i; j < N; j += i) { if (mu[j / i] == 1) g[j] = ((ll)g[j] * f[i]) % MOD; else if (mu[j / i] == -1) g[j] = ((ll)g[j] * qmi(f[i], MOD - 2)) % MOD; } } for (int i = 1; i < N; i++) g[i] = (ll)g[i] * g[i - 1] % MOD; }
int main() { init(); int T; scanf("%d", &T); while (T--) { scanf("%d%d", &n, &m); int k = min(n, m); int ans = 1; for (int l = 1, r; l <= k; l = r + 1) { r = min(n / (n / l), m / (m / l)); int t = (ll)(n / l) * (m / l) % (MOD - 1); ans = ((ll)ans * qmi((ll)g[r] * qmi(g[l - 1], MOD - 2) % MOD, t)) % MOD; } printf("%d\n", ((ll)ans + MOD) % MOD); } return 0; }
|